#This script will determine the residual value of the relationship between ES and a species abundance (tau=.975)
#Clear workspace

rm(list=ls())

in.dir = "C:/R/In/"
out.dir = "C:/R/Out/"

#Set working directory

setwd(in.dir)                 

#load necessary libraries

library("SDMTools")
library("quantreg")
library("adehabitat")

#Import the dataset for analysis
master=NULL
master=samples= read.csv("samplesmaster.csv",header=T)
s.fields = names(samples)

#Perform quantile regression for each of 4 species against it's own ES values

cr.rq=rq(samples$cr_abund~samples$cr_es,data=samples,tau=.975,method="br",model=T)
sb.rq=rq(samples$sb_abund~samples$sb_es,data=samples,tau=.975,method="br",model=T)
gq.rq=rq(samples$gq_abund~samples$gq_es,data=samples,tau=.975,method="br",model=T)
lc.rq=rq(samples$lc_abund~samples$lc_es,data=samples,tau=.975,method="br",model=T)

#Plot phys data

png(paste(out.dir,"sbcompnew.png", sep=""), height=1000, width=1000)
par(cex=c(1.5),cex.lab=c(1.5))
#plot(master$airtemp,master$sb.rq.residuals,xlab="AIR TEMP", ylab="ABUNDANCE DIFFERENCE", main="SAPBASI",ylim=c(-5,7.5))
#plot(master$airtemp,master$lc.rq.residuals,xlab="AIR TEMP", ylab="ABUNDANCE DIFFERENCE", main="LAMCOGG")
#plot(master$airtemp,master$gq.rq.residuals,xlab="AIR TEMP", ylab="ABUNDANCE DIFFERENCE", main="Prickly Forest Skink")
#plot(master$hum,master$gq.rq.residuals,xlab="HUMIDITY", ylab="ABUNDANCE DIFFERENCE", main="PRICKLY FOREST SKINK")
plot(master$lc_abund, master$sb.rq.residuals, xlab="ABUNDANCE OF COMPETITOR", ylab="ABUNDANCE DIFFERENCE", main="SAPROSCINCUS BASILISCUS", ylim=c(-5,4), xlim=c(0,70))

dev.off()

#Make a data frame from the regression residuals and the raw georef data

master=data.frame(samples, cr.rq$residuals, gq.rq$residuals,lc.rq$residuals, sb.rq$residuals)
samples=master
s.fields = names(samples)

 #First begin by working with the samples data, plotting relationships between all variables and sample abundance
#Use a for loop to perform both quantile and OLS regression between all variables and samples abundance

num=1

for(fields in s.fields)

  {
  
  t.data=data.frame(samples$cr.rq.residuals, samples[,num])

  png(paste(out.dir,"CARRUBR","/",fields,".png", sep=""), height=800, width=800)

  plot(t.data[,2], t.data[,1], xlab=fields, ylab="residuals", main="CARRUBR")
  
  t.rq1 = rq(t.data[,1]~t.data[,2], tau=.90, data=t.data, method="br", model=T)
  t.rq2 = rq(t.data[,1]~t.data[,2], tau=.95, data=t.data, method="br", model=T)
  t.rq3 = rq(t.data[,1]~t.data[,2], tau=.975, data=t.data, method="br", model=T)
  t.rq4 = rq(t.data[,1]~t.data[,2], tau=.99, data=t.data, method="br", model=T)
  t.ols = lm(t.data[,1]~t.data[,2], data=t.data)

  abline(t.rq1, col="red")
  abline(t.rq2, col="red")
  abline(t.rq3, col="red")
  abline(t.rq4, col="red")
  abline(t.ols, col="blue")

  
  dev.off()
  
  num=num+1

  }
  
  num=1
  for(fields in s.fields)

  {

  t.data=data.frame(samples$sb.rq.residuals, samples[,num])

  png(paste(out.dir,"SAPBASI","/",fields,".png", sep=""), height=800, width=800)

  plot(t.data[,2], t.data[,1], xlab=fields, ylab="residual", main="SAPBASI")

  t.rq1 = rq(t.data[,1]~t.data[,2], tau=.90, data=t.data, method="br", model=T)
  t.rq2 = rq(t.data[,1]~t.data[,2], tau=.95, data=t.data, method="br", model=T)
  t.rq3 = rq(t.data[,1]~t.data[,2], tau=.975, data=t.data, method="br", model=T)
  t.rq4 = rq(t.data[,1]~t.data[,2], tau=.99, data=t.data, method="br", model=T)
  t.ols = lm(t.data[,1]~t.data[,2], data=t.data)

  abline(t.rq1, col="red")
  abline(t.rq2, col="red")
  abline(t.rq3, col="red")
  abline(t.rq4, col="red")
  abline(t.ols, col="blue")


  dev.off()

  num=num+1

  }

  num=1
  for(fields in s.fields)

  {

  t.data=data.frame(samples$gq.rq.residuals, samples[,num])

  png(paste(out.dir,"GNYQUEE","/",fields,".png", sep=""), height=800, width=800)

  plot(t.data[,2], t.data[,1], xlab=fields, ylab="residuals", main="GNYQUEE")

  t.rq1 = rq(t.data[,1]~t.data[,2], tau=.90, data=t.data, method="br", model=T)
  t.rq2 = rq(t.data[,1]~t.data[,2], tau=.95, data=t.data, method="br", model=T)
  t.rq3 = rq(t.data[,1]~t.data[,2], tau=.975, data=t.data, method="br", model=T)
  t.rq4 = rq(t.data[,1]~t.data[,2], tau=.99, data=t.data, method="br", model=T)
  t.ols = lm(t.data[,1]~t.data[,2], data=t.data)

  abline(t.rq1, col="red")
  abline(t.rq2, col="red")
  abline(t.rq3, col="red")
  abline(t.rq4, col="red")
  abline(t.ols, col="blue")


  dev.off()

  num=num+1

  }
  
  num=1
  for(fields in s.fields)

  {

  t.data=data.frame(samples$lc.rq.residuals, samples[,num])

  png(paste(out.dir,"LAMCOGG","/",fields,".png", sep=""), height=800, width=800)

  plot(t.data[,2], t.data[,1], xlab=fields, ylab="residuals", main="LAMCOGG")

  t.rq1 = rq(t.data[,1]~t.data[,2], tau=.90, data=t.data, method="br", model=T)
  t.rq2 = rq(t.data[,1]~t.data[,2], tau=.95, data=t.data, method="br", model=T)
  t.rq3 = rq(t.data[,1]~t.data[,2], tau=.975, data=t.data, method="br", model=T)
  t.rq4 = rq(t.data[,1]~t.data[,2], tau=.99, data=t.data, method="br", model=T)
  t.ols = lm(t.data[,1]~t.data[,2], data=t.data)

  abline(t.rq1, col="red")
  abline(t.rq2, col="red")
  abline(t.rq3, col="red")
  abline(t.rq4, col="red")
  abline(t.ols, col="blue")


  dev.off()

  num=num+1

  }
  
#Write out this data frame as a .csv file

write.csv(x=master, file=paste(out.dir,"quantregresidsbygeoref.csv",sep""),row.names=F)

#Perform a multiple step wise regression of residual values against HEMI, coastdist, dem80, and slope

cr.lm = lm(master$cr.rq.residuals~master$canopen_road+master$coastdist+master$dem80+master$LAI4_forest+master$LAI4_road+master$MODIS_mean+master$MODIS_sd+master$slope, data=master)
cr2.lm = lm(master$cr.rq.residuals~master$canopen_road+master$dem80+master$LAI4_road+master$bc05micro+master$bcdiff, data=master)
cr3.lm = lm(master$cr.rq.residuals~master$bc05micro+master$bcdiff, data=master)

gq.lm = lm(master$gq.rq.residuals~master$canopen_road+master$coastdist+master$dem80+master$LAI4_forest+master$LAI4_road+master$MODIS_mean+master$MODIS_sd+master$slope, data=master)
gq2.lm = lm(master$gq.rq.residuals~master$canopen_road+master$dem80+master$LAI4_road+master$bc05micro+master$bcdiff, data=master)
gq3.lm = lm(master$gq.rq.residuals~master$bc05micro+master$bcdiff, data=master)

lc.lm = lm(master$lc.rq.residuals~master$lsc_mean+master$canopen_road+master$coastdist+master$dem80+master$LAI4_forest+master$LAI4_road+master$MODIS_mean+master$MODIS_sd+master$slope, data=master)
lc2.lm = lm(master$lc.rq.residuals~master$lsc_mean+master$canopen_road+master$dem80+master$LAI4_road+master$bc05micro+master$bcdiff, data=master)
lc3.lm = lm(master$lc.rq.residuals~master$bc05micro+master$bcdiff, data=master)

sb.lm = lm(master$sb.rq.residuals~master$lsc_mean+master$canopen_road+master$coastdist+master$dem80+master$LAI4_forest+master$LAI4_road+master$MODIS_mean+master$MODIS_sd+master$slope, data=master)
sb2.lm = lm(master$sb.rq.residuals~master$lsc_mean+master$canopen_road+master$dem80+master$LAI4_road+master$bc05micro+master$bcdiff, data=master)
sb3.lm = lm(master$sb.rq.residuals~master$bc05micro+master$bcdiff, data=master)
